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Abstract 

The matter-enhanced oscillations of two neutrino flavors are studied using 
a uniform semiclassical approximation. Unlike some analytic studies which 
have focused on certain exactly-solvable densities, this method can be used 
for an arbitrary monotonic density profile. The method is applicable to a 
wider range of mixing parameters than previous approximate methods for 
arbitrary densities. The approximation is excellent in the adiabatic regime 
and up to the extreme nonadiabatic limit. In particular, the range of validity 
for this approximation extends farther into the nonadiabatic regime than for 
the linear Landau-Zener result. This method also allows calculation of the 
source- and detector-dependent terms in the unaveraged survival probability, 
and analytic results for these terms are given. These interference terms may 
be important in studying neutrino mixing in the sun or in supernovae. 
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I. INTRODUCTION 



Matter-enhanced oscillations of neutrino flavors via the Mikheyev-Smirnov-Wolfenstein 
(MSW) mechanism M have been studied for neutrinos in various environments, but most 
extensively for the sun, in connection with the solar neutrino problem 0. For a recent 
review of the solar neutrino problem and the ongoing neutrino detection experiments, see 
Ref. [[J. Recently, interest has also been developing for the study of neutrino oscillations in 
super novae |4]]. 

The approximate results derived in this paper are applicable to matter-enhanced two- 
flavor neutrino oscillations in general physical situations. Analytic results are important for 
several reasons. While numerical integration of the MSW equations is straightforward, it 
becomes extremely tedious when it must be done for a large range of the mixing parameters. 
Analytic results also allow a greater understanding of the effects of changes in the parameters, 
and may be useful for extracting information about the solar density from the measured 
neutrino fluxes. 

Analytic studies of matter-enhanced neutrino oscillations proceed along two lines. The 
first approach is the study of certain densities for which an exact solution for the oscillation 
probability can be obtained. The mixing parameters are allowed to be arbitrary. The 
exponential density has attracted particular interest, since it approximates the solar density. 
A catalog of all of the exactly solvable densities has been presented in Ref. |J. The second 
approach allows for a general density, but restricts the parameters so that an approximation 
can be made to the equations of motion, which are then solved exactly. The approximations 
are chosen so that the exact results are recovered in either the extreme nonadiabatic or 
extreme adiabatic limit. In this paper, we consider a uniform semi classical approximation 
to derive the neutrino conversion probability for an arbitrary density. The solution is exact 
in the adiabatic limit, like the linear Landau-Zener result. However, the new result has a 
larger range of validity in the nonadiabatic regime. In the body of the paper, we will discuss 
how some of the different approximations are related. 



II. MATTER-ENHANCED NEUTRINO OSCILLATIONS 

A. Coupled Equations in the Flavor Basis 

For two neutrino flavors (taken here to be electron and muon) in matter, the equations 
of motion for the v e and probability amplitudes in the relativistic limit are 
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A — 5m 2 cos 29 v Sm 2 sin 29 v 
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5m 2 sin 29 v —A + 5m 2 cos 9 V 







(2.1) 



where all terms in the Hamiltonian proportional to the identity have been dropped since 
they do not contribute to the relative phase between the v e and z/ M components. The vacuum 
mixing parameters are specified by the vacuum mixing angle 9 V , taken to be < 9 V < 7r/4, 
and the vacuum mass-squared splitting 5m 2 = m\ — m 2 , where we take m 2 > mi. Electron 
neutrinos experience charged- current scattering with the electrons in the medium, whereas 
muon neutrinos do not. This difference yields the effective mass correction 
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A = 2V2 G F N e (t)E , 



(2.2) 



where Gf is the Fermi constant and N e (t) is the number density of electrons in the medium. 

Before proceeding further, we switch to working with dimensionless quantities. We define 
a length scale 



5m? /IE 



(2.3) 



and use this to define x = t/L. Since we will be making a semiclassical expansion, we need 
to be able to keep track of formal powers of H. For each h in the problem, we write A 
and consider A to be formally small; this is equivalent to saying that the length L is small. 
We will make expansions in powers of A, truncating the higher orders. At the end of the 
calculation, we will set A = 1. For notational convenience, we write 



zA— 

dx 



We have defined 



" * e (x) " 
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= H e ^(x) 















r](p(x) vA 

A —7]ip(x) 



(2.4) 



f]ip(x) = ((x) — co$26 v 



and 



A = sin 26 v . 

The scaled electron density is 

C(x) = C(x,)A e (x)/A e (x,) 
normalized at the initial point 

2V2 GfENJxA 



5m 2 



(2.5) 



(2.6) 



(2.7) 



(2.8) 



Note that there are notation changes from Ref. |§-§|; here we have made A and if dimen- 
sionless. The factor 77 (taken to be ±1), is introduced above to control the analytic behavior 
of the function <f(x) in the complex plane, as explained in the Appendix. In the expressions 
with ip 2 below, we drop rf = 1. 



B. Coupled Equations in the Adiabatic Basis 



The flavor-basis Hamiltonian of Eq. ( j2.4| ) can be instantaneously diagonalized. We make 
a change of basis 



#i(a;) " 




' tf e (a:) " 






= R(-6(x)) 






* 2 (x) _ 









cos6(x) —sm6(x) 
sin^(x) cos^(x) 



(2.9) 
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\E'i(x) is the probability amplitude to be in the "light" (primarily electron type) eigenstate in 
the mass basis, and \E r 2 (a : ) is the probability amplitude to be in the "heavy" (primarily muon 
type) eigenstate in the mass basis. The requirement that this transformation instantaneously 
diagonalize H efJL (x) defines the matter mixing angle via 



sm29(x) 



and 



cos 29(x) 



A + (p 2 (x) 



(2.10) 



(2.11) 



The matter angle thus ranges from tt/2 at infinite density to 9 V in vacuum. At the resonance, 
9 = 7r/4. The instantaneous eigenvalues of H e/I (x) are 



=F JA + <p 2 (x) 



(2.12) 



corresponding to (the "light" eigenstate) and ^2 (the "heavy" eigenstate), respectively. 
The splitting between the instantaneous mass eigenstates has a minimum as a function of x 
when <^(x) = 0, or ((x) = cos 29 v \ this is the MSW resonance point, which will be denoted by 
x c . The trajectories of these eigenvalues represent an avoided level crossing. The adiabatic 
limit is the case where the neutrino stays in one of the instantaneous eigenstates during its 
entire propagation. In the nonadiabatic limit, the neutrino may "hop" from one eigenstate 
to the other near the resonance. 

In the mass basis, the equations of motion are 



iX 



d_ 

dx 



'*i(x)~ 










= H 12 {x) 






V 2 {x) _ 









-^/A + ip 2 (x) -i\9'(x) 
i\9'{x) ^A + <p 2 (x) 



[2.13) 



Throughout the paper, prime denotes derivative with respect to x. When the density is 
changing slowly, then so is the matter angle 9(x), and the off-diagonal terms can be neglected; 
for that reason, this is also known as the "adiabatic" basis. The adiabaticity parameter is 
defined as 



70) 



V / A + ^ 2 (^) 



iX9'(x) 



(2.14) 



where 9'(x) can be derived from Eqs. ( |2.10| ) and (|2 . 1 1| ) . When this parameter ^{x) is large, 
we can neglect the off-diagonal terms. All nonadiabatic behavior, i.e., hopping from one 
mass eigenstate to the other, takes place in a neighborhood of the resonance. It is there that 
~f(x) is minimized, so the requirement of 7(x) ^> 1 for adiabatic propagation is the most 
exacting: 



7c = l{x c ) 



2 sin 2 26, 



1 



A cos 29 v IC'/CU 



> 1. 



(2.15) 
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In this limit, the equations of motion can be integrated immediately, yielding pure phases 
for ^i(x) and ^(e). At the initial point X{, we take the neutrino to be a pure u e , so 



tte(Xi) = 1 (2.16) 



and 



dVJx) 



dx 



\rypixi). (2.17) 



The latter follows from *ff e (xi) = 1, ^(x;) = 0. We denote the initial matter angle by 9i 
and reference the phases from the initial point Xi. The phase integral will be denoted by 

I p (x, Xi) = [ X dxJk + (p 2 (x) . (2.18) 

JXi 

Taking into account the basis changes at the initial and final points, the adiabatic solutions, 
valid in the limit "f(x c ) ^> 1, are 

\l/ e (x) = cos 6{x) cos 6i exp ^+—I p (x,Xi)^ + sin 6{x) sin 6>j exp ^— — I p (x, Xi)^j (2.19) 

and 

^n( x ) — — srn cos 0i ex P (+~^Ip( x > x i)^j + cos Q( x ) srn Qi ex P -^Ip( x i x i)^J ■ (2.20) 

These forms hold both before and after the resonance in the adiabatic limit. If nonadiabatic 
corrections are taken into account, then the wave functions will have these forms before the 
resonance but will be more complicated after the resonance. In the adiabatic limit ||, the 
electron neutrino survival probability at a general point x is: 

P{y e -> u e )(x,Xi) = \^ e ( x )\ 2 

= - [1 + cos 29i cos 20(e)] + r sin 29 i sin 29(x) cos ( -I p {x, Xi)) . (2.21) 
2 2 \ A / 

Note that the second term depends upon the source and detector positions, and will dis- 
appear under averaging of either. The probability of conversion to muon type is given by 
P[y e — > Vy) = 1 — P[y e — > u e ). If the final point is chosen in vacuum, then 9(x) — > 9 V . 

With the adiabatic limit in hand, the obvious thing to do is to seek the corrections that 
take into account Phop, the probability of hopping from one mass eigenstate to the other. 
Above, the adiabatic approximation was controlled by the ratio of diagonal to off-diagonal 
elements. That ratio is in turn controlled by A, which keeps track of powers of h. In the 
semiclassical limit of h — > 0, one has A — > and 7 C — > oo. Note that A appears above in the 
adiabatic survival probability only in the phase; the fully averaged expression is independent 
of A. The way to treat Phop systematically is to expand in powers of A and to keep only the 
lowest-order terms. 
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C. Uncoupled Equations in the Flavor Basis 



The coupled first-order equations for the flavor-basis wave functions can be decoupled to 
yield 



A 



X 



dx 2 



A + <p 2 (x) + i\r](p'(x)} * e (x) = 



and 



A 



2 d 2 ^,(x) 



dx" 1 



A + p 2 (x) — i\r)ip'(x) 







(2.22) 



(2.23) 



2~5| ) and (|2.6|). Such a simple decoupling is not 



where (f(x) and A are defined in Eqs. 
possible in the matter basis. 

These Schrodinger-like equations are similar to those for non-relativistic particles in the 
presence of a complex barrier, and for convenience we use the language of wave mechanics 
to describe them. In particular, to the extent that we can ignore the imaginary terms in 
the potential, these correspond to particles above a barrier (since A > 0). There are two 
caveats regarding discussing this as a barrier penetration problem. First, that our boundary 
conditions do not correspond to the usual picture of incident, reflected, and transmitted 
waves; in general, there are waves moving in each direction on each side of the barrier. 
Second, the pure imaginary terms in the potentials play an extremely important role here, 
even in the asymptotic regions. These terms are needed not only to represent nonadiabatic 
transitions, but also to keep up with the local matter angle. 

In this problem, then, the quantity of interest is not a reflection or transmission coeffi- 
cient, but rather P[y e — > v e ) = |\l/ e (x — > oo)| 2 , the probability of the neutrino being of the 
electron type far from the source. In general, this is a function of both source and detector 
positions, though typically only the fully averaged result is presented. However, those inter- 
ference terms could be important, and we present approximate expressions for them in the 
next section. 

Two well-known semiclassical treatments of this problem are via the Wentzel-Kramers- 
Brillouin ( WKB) and linear Landau-Zener [10-12H methods. The WKB technique globally 
maps the "potential" discussed above onto the free-particle potential (i.e., a constant den- 
sity). By "global mapping", we mean a variable stretching of the axis that deforms the shape 
of one potential into another. In fact, the WKB treatment turns out to be identical to the 
adiabatic approximation 0. The linear Landau-Zener technique locally maps onto a linear 
density (i.e., extends a linear profile from a single -MSW resonance- point with the right 
density and derivative), and hence a "potential" with a parabolic real part and constant 
imaginary part. While the linear Landau-Zener result is easy to derive and apply, there 
are two problems. First, it is notoriously difficult to get the boundary conditions right (for 
a complete explanation of how to handle this, see Ref. Hl3| ). Second, since Landau-Zener 
is a point mapping, the expression for Phop is not very accurate. In the case of neutrino 
oscillations in of the sun, the exponential Landau-Zener approximation circumvents these 



problems | I3fl . 

The aim of this paper is to calculate the nonadiabatic corrections semiclassically, but 
with a global mapping of the "potential", where as in the Landau-Zener calculation we 
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choose as a model the case of a linear density. By using a global mapping, the correct 
boundary conditions are automatic. Further, the expression for Ph op is more accurate. The 
approximate wave function is uniformly valid in x (though the approximation is not uniform 
in the mixing parameters). 



III. UNIFORM SEMICLASSICAL SOLUTION OF THE MSW EQUATIONS 

A. Semiclassical Background 

In the adiabatic limit, only the lowest order is kept in the limit A — > in Eq. ( |2.13| ), so 
the Hamiltonian is taken as diagonal (no hopping from one mass eigenstate to the other) 
and the integration is trivial. The treatment at that order suggests that in order to take 
into account the probability of hopping, we will need to consider further orders in A. In this 
section, we will show that it is possible to obtain a rather accurate expression for the electron 
neutrino survival probability by making a semiclassical expansion, i.e., by considering only 
the two lowest orders in A when solving the MSW equations. 

The expressions derived below will hold for values of the mixing parameters from the 
extreme adiabatic limit up until the extreme nonadiabatic limit. In order to obtain solutions 
that hold in the extreme nonadiabatic limit, one would formally have to consider all orders 
in A. Since semiclassical expansions are asymptotic (i.e., non-convergent) in general, it is not 
clear that this would work in practice. A much better approach for the extreme nonadiabatic 
limit is to consider expansions in 1 / A |14[ . 



Semiclassical methods (for reviews, see Ref. |15|) are used in quantum mechanics to 
provide approximate solutions to the Schrodinger equation in the limit that A is small. As 
noted, in the WKB method, one bases the approximate solutions on free-particle solutions. 



A procedure was developed by Miller and Good Jl6|] that instead bases the approximate 
solutions on the known solutions of a Schrodinger equation with a similar potential. In this 
method, the turning point singularities of the primitive WKB method are regulated, and 
the solutions are uniformly valid: they hold over the whole range in x and are well-behaved 
at the turning points. A further advantage of the Miller-Good method is that different 
potentials are treated with the same formalism, i.e., the method of connection is the same 
for all potentials with the same number and type of turning points. 

The MSW equations of Eqs. ( |2.22|) and ( j2.23| ) are Schrodinger-type equations for particles 



in the presence of complex potentials of the form V(x) = — [f 2 {x) ± i\f]ip'(x)}, with (p(x) 
independent of A. In the Appendix, we summarize the extension of the uniform semiclassical 
approximation to treat potentials with this specific dependence on A, originally introduced in 
Ref. [0-0 • This special form of the potential arises in supersymmetric quantum mechanics; 
see Ref. IHI for discussion. 



B. Application to the MSW Problem 

The method presented in the Appendix allows a uniform semiclassical solution for \l/ e (x). 
By "uniform", we mean that the local error incurred by using the approximate solution 
developed there in the exact differential equation is bounded as a function of x. This is to be 
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distinguished from a semiclassical solution via the primitive WKB method. Such as solution 
has an unbounded error near a turning point. For MSW propagation, the turning points are 
complex (they are near the resonance point). In the nonadiabatic limit, the turning points 
approach the real axis, which means that the WKB solutions are unable to represent any of 
the nonadiabatic behavior. In contrast, the uniform semiclassical approximation used here 
is excellent for all but the extreme nonadiabatic limit. Since we will make a semiclassical 
expansion, we explicitly show all factors of h (via A). Either an increasing or decreasing 
density can be considered, by proper choice of r\. 

From the derivation in the Appendix, the general solution to Eq. ( |2.22| ) is 

* e (x) = K(x) [AD v {z{x)) + BD v {-z{x))\ , (3.1) 

where 

y= "- 1 - gl /\ (3.2) 

with 

Q = - p dxJA + if 2 (x) . (3.3) 

The limits of integration above are the zeroes of the integrand, chosen as described in the 
Appendix. The argument of the parabolic cylinder functions is given by 

z(x) = S(x) » ^ (S (x) + XS 1 (x)) , (3.4) 

where Sq(x) and Si(x) are described in the Appendix. 

Given appropriate initial conditions, one can solve for A and B. In some situations, it 
may be useful to evaluate \I/ e (a;) for all x. This requires evaluating the gamma function 
for complex argument and the parabolic cylinder functions for general complex order and 
argument. For general comments on routines available for the numerical evaluation of special 
functions, see Ref. ( [[L7[|). The gamma function for complex argument can be evaluated 
with CERNLIB fI8| . General properties of the parabolic cylinder functions may be found 
in Ref. [p!9| -|22]|. While library routines do exist for various special cases of the parabolic 



cylinder functions, to our knowledge there is nothing available that is general enough p3| 



The technique^ used here is to use the power series |22| for small \z\, the asymptotic series 



pOfl for large \z\, and direct numerical integration of the defining differential equation with 
ODEPACK [24] for moderate Fortunately, one does not generally have to perform any 
integrations for the parabolic cylinder functions, as only the asymptotic forms are needed. 

We will use the asymptotic forms at both the production and detection points. As shown 
below, this means that we assume adiabatic propagation at those two points. This matching 
is justified to the extent that the production and detection points are sufficiently far from the 
resonance point. In practice, these requirements do not present any difficulties. Consider 



x The code is available upon request from the authors. 
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the sun as an example, with neutrinos produced at the solar center. If the resonance is 
near the production point, or there is no resonance, then this implies that 5m? is large 
and the entire propagation is adiabatic, except for extremely small mixing angles. If the 
resonance is at a very low density, i.e., approaching vacuum, then 8m? is very small and 
our approximation breaks down for other reasons described below. Note that in the linear 
Landau-Zener treatment, one has to handle the final x point carefully since the density runs 
negative at large x, and 9(x) — > 0, not 9 V . No such difficulties with the boundary value of 
the matter angle arise in the treatment given here. 

The asymptotic forms developed below represent ^ e (x) well for large but finite x. All 
of the expansions below are just to get outside of the resonance region; we do not take x 
so large that the matter angle is either 7r/2 or 9 V . More precisely, (see the discussion in the 
Appendix), the approximate solutions are characterized by two scales, one set by Sq(x) and 
the other by (p(x). The function S (x) is asymptotic outside the resonance region, whereas 
<p(x) is not asymptotic until the density is zero or infinite. In this formulation, S (x), but 
not tp(x), will be taken to be asymptotic. This means that we have the control to connect 
opposite sides of the resonance region without having to take x —>■ ±oo, i.e., we do not have 
to extend the density profile indefinitely. 



C. Asymptotic Solutions and Connection Formulae 

Using the definition of the matter angle given in Eqs. ( [2.10| ) and ( |2.11| ), we can rewrite the 
pre-exponential factors in the asymptotic solution of \I/ e (x), Eq. (|A37|) . There are two cases, 
depending on how (p and r] are chosen. The first case has tp(x) = cos 29 v — ((x),r] = —1, so 



A 



A(A + ip 2 (x)) 



1/4 



if + VA + if 2 



-1/2 



cos9(x) , 



(3.5) 



and 



A 



nl/4 



4(A + <p 2 (x)) 



> + VAT?' 



"1/2 



A 



sin9(x) . 



(3.6) 



In the other case of tp(x) = cos 29 v — ((x),r] = —1, these are reversed. In either case, the 
prefactors associated with the various terms are: 

C_ : cos 6* (a;) 
D- : sin#(x) 
C + : sin 9 (x) 



cosy [x\ 



(3.7) 



It is important to note that these will be evaluated at general values of x outside the 
resonance region; |x| will not be so large that 9{x) — > n/2 or 9(x) — * 9 V . The asymptotic 
wave functions are 



\l/ e (x — > — oo) = C_ cos6'(a;) exp ( +— I p (x, Xi) ) + -D_ sin#(a;) exp ( ——I p (x, x 



A 



(3.8) 
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and 

^ e (x — > +00) = C + sin 9{x) exp ^— — J p (x, + D + cos6>(x) exp (+—I p (x,Xi)^ . (3.9) 

From their form, we see immediately that the asymptotic wave functions represent adiabatic 
propagation. The coefficients C±,D± still depend on r], which will allow us to consider 
increasing or decreasing densities. The phase integral function I p is defined in Eq. QA12 ). 

With the asymptotic wave function in this form, it is rather easy to apply the initial 
conditions. As before, we take the neutrino at the initial point Xi to be a pure u e , so 

*e(Xi) = 1 (3.10) 

and 

i 



Tl . . dV K (x) 



dx 



jWixi). (3.11) 



In either case regarding the signs of ip and 77, one immediately obtains C_ = cos ft and 
D- = sin ft, so 

ty e (x — > —00) = cos9(x) cos ft; exp (+^I p (x, Xi)j + sin 6 (x) sin ft exp ^— — I p (x, Xi) \ . 

(3.12) 

which is of course the adiabatic solution given in Sec. ( |I B| ). 

We now turn to the evaluation of the coefficients C + and D + that are needed after the 
resonance. From the above and Eqs. (|A39| ) and ( |A41|) , 

C_ = cos ft = Cexp (+ l -ReI p (x u x Q )^j (Ae~ wn + B) (3.13) 

D_ = sin ft = Cexp (-jReI p (x h x )) (Ae~ W7T ) . (3.14) 
These determine the coefficients A and B of the general solution: 

A = IT 1 exp (+jReI p {xi, x )) e^sinft (3.15) 

B = C^ 1 exp ^— — Kel p (xi, x )^ cos ft — D^ 1 exp Re/ p (xj, x )^ sin ft. (3.16) 

Then 

C + = 2iCD~ 1 e iU7T sin 6i + exp (-jReI p {xi, x f) e"^ cos ft (3.17) 

D + = C-^e-^cosft - exp (+jReI p (xi, x )) e^" sin ft , (3.18) 

where C and .D are given in Eqs. ( |A42j) and ( |A43| ) . 

The asymptotic forms of ^ e (x) shown above are perfectly general, and depend only on 
the assumption of adiabatic propagation outside the resonance region. The heart of this 
problem is the connection of the asymptotic coefficients C_ and D- to C + and D + . That 
connection represents the integration of the solutions through the resonance region. In our 
case, that information is carried by the coefficients A and B of the general (but approximate, 
due to the mapping) solution in terms of parabolic cylinder functions. 
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D. Resonance Transition Coefficients 



Above, the asymptotic wave functions were written in terms of the adiabatic solutions, 
which is convenient for applying initial conditions and deducing the connection formulae. 
Before squaring the asymptotic wave function to obtain the neutrino survival probability, it 
is convenient to rewrite the wave function in a slightly different form: 



+00) = C + sin 9{x) exp ^— — I p (x, Xi)^J + D + cos9(x) exp \+—I p (x,Xi)\ 
C\ sin 9i exp (+—ReI p (xi, x )^ + c 2 cos 6^ exp ^— — Rel p (xi, x ] 



x sin 9(x) exp ^— — Rel p (x, x )^j 



+ 



cos 6^ exp ( — — Rel p (xi, 20)^) + d 2 sin 6^ exp (+—ReI p (xi, xq] 



x cos 6{x) exp ( +— Rel p (x, xq\ 



(3.19) 



The terms inside the square brackets depend only on the source position 2j, whereas the 
terms outside depend only on the final position x. Since all of the adiabatic phases and 
matter angles for the asymptotic solutions are written explicitly, the matrix of coefficients 
given by Ci, C2, g?i, (£2 represents only the nonadiabatic transitions in the resonance region. 
These coefficients change in the resonance region, but tend to asymptotically to constants 
outside of it. Since the 2x2 Hamiltonian is Hermitian and traceless, the time-evolution 
operator must be a member of the SU(2) group. Thus this matrix must assume the form 



Cl c 2 



-Co 



(3.20) 



where |ci| 2 + |c 2 | 2 = 1. 

By comparison to the forms of C + and D + in Eqs. QA38 ) and ( |A40| ) , the new coefficients 
are easily found to be 



Cl 



c 2 = e 



2vr VA 



( iQ\ . , 



(3.21) 
(3.22) 



2^F /fi X+^/2A-,/2 / e -«r/2\-" 



r(-iz) V a 



d 2 = — e" 



,+3mt/ 



/ 4 V2exp ( ~ 



(3.23) 
(3.24) 



By analysis of two cases of 77 = ±1 separately, one can show 



r -1/ r — 



7T 



V2A/ sinh (Qtt/2X) ' 



(3.25) 
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2 T(-v)\ 2 ( | sin(z/7r)| V^ 2A = 1 - e~^' 2X , (3.26) 



7T 



where v is given by Eq. ( |A22|) . Using these relations, it may easily be verified explicitly that 
d\ = c\ and d* 2 = — c 2 , and that 



\ c i\ 



1 - e~ Q7T/x , (3.27) 



\ C 2\ 



,-flTT/X 



(3.28) 



so |ci| 2 + |c 2 | 2 



1. 



Starting with Eq. ( |2.13p , one can determine how C\ and c 2 depend on 77, i.e., on whether 
the density is increasing or decreasing. One can show that C\ must be independent of 77, and 
that c 2 must change sign if 77 does. With the present form of c\, this is not obvious. Define 
the phase of c\ as 



ci 



\ci\e 



(3.29) 



An asymptotic series can be developed for this phase a. Using the special form of the Stirling 
expansion of the gamma function for purely imaginary argument given in Eqs. (6.1.43-44) 



of Ref . |l9|] , one can show that the phase of c\ , in the limit Vl/ A is large, is 



a 



E 



:-i) n - i B 2n (2\ 



2n-l 



{ 2n(2n - I) 



(3.30) 



where B 2n are Bernoulli numbers 19[| . When a linear density is considered, this expression 
for the phase is equal to that given in Ref. [|25[| . We do not use this limit for the phase 
in general, since it requires that Q/X > 1, which is unnecessarily restrictive on the range 
of validity of our main approximation. Since this is independent of 77, so is c\. Note from 
the definition of v in Eq. ( |A22| ) that c 2 = e~ W7r is in fact a real number, though it may be 
positive or negative, and changes sign if r\ does. 



E. Calculation of P{y e — > v e 



The electron neutrino survival probability at a general point x after the resonance for a 
neutrino produced at Xj is given by the modulus squared of the amplitude \l/ e (x — > +00) for 
the neutrino to be of the electron type. First write 



-00 



Ci sin 6i exp y+—ReI p (x i} x ] 
x sin 9(x) exp ^— — Re/ P (a;, x )^j 



+ c* cosOi exp y--ReI p (x il x ] 
x cos 9(x) exp ^+—Relp(x,x )^j 



c 2 cos 6i exp ( --Rel p (xi, x \ 



c* 2 sin Oi exp ( +-ReI p (x il x ] 



(3.31) 
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After taking the squared modulus of this expression for ty e (x), and then reducing it, the 
survival probability takes the form 



P(u e -> u e )(x,Xi) 



- [l+ (1 -2|c 2 | 2 )cos2^cos2^(x 

1 , (2 
'c 1 |c 2 sin 29 i cos26>(a;) cos ( —Rel p (xi, x ) + a 



A 



+ - sin 29 i sin 2#(x) cos ( — Rel p (x, xi) — 2a 
2 V A 

— |c 2 | 2 sin26 l jsin26'(a;) cos ^— Rel p (xi, x ) + «^ cos ^— Re/ p (a;, x ) — 



a 



+ | Ci | c 2 cos 29 i sin 26 , (x) cos ( — Rel p (x, x ) — a 



(3.32) 



This, along with Eqs. fl3.21|) and (|3.22 ), is our main result. Recall that in our approximation, 
c 2 is real. In general, the phase a should be extracted from c± directly, rather than taken from 
the asymptotic series for a given above. This simple form for the survival probability can be 
evaluated easily and rapidly, providing accurate results for both the direct and interference 
terms for all mixing parameters except for the extreme nonadiabatic limit. It is much more 
convenient than direct numerical solution of the MSW equations, especially if many values 
of the mixing parameters need to be explored. 



When n/X 3> 1, i.e., the adiabatic limit, \ci\ 
form for the survival probability reduces to 



1, |c 2 | — > 0, a — > 0, and this general 



P(u e — > v e )(x, Xi) —> — [! + cos 29{ cos 29{x)\ H — sin 26^ sin2^(x) cos ( — Rel p (x, x 



A 



(3.33) 



which is the usual adiabatic result. 



Typically, the final point x will be taken in vacuum, so 9{x) — > 9 V and JA + cp 2 (x 



and 



2i 

exp ( ±— Re/ P (a;, xq] 



const, x exp (±2ix/\) . 



(3.34) 



(The same applies when the lower limit in I p is x«, though the constant will be different.) 
The oscillation length in vacuum is ttL, where L is given by Eq. (2T3). For example, in 
the solar neutrino problem, the favored MSW parameters lead to an oscillation length of 
~ 1000 km p6 |. In such cases, where the oscillation length in vacuum is much less than the 
variation in the source-detector distance, it will be appropriate to average over the detector 
position. If that is done, then the survival probability no longer depends on x, but does still 
depend on X{ and is given by 



P(v P 



u e )(xi) = - [l + (1 - 2|c 2 | 2 ) cos 29 { cos 29 v 

— — | ci | C2 sin 29i cos 29 v cos ^— Rel p (xi, Xq) + a ) . 



(3.35) 



This shows that the source term may be important even after detector averaging. If the 
source is extended, or if an energy spectrum is considered, one can also average over the 
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source position. The completely averaged result for the electron survival probability is then 
given by 



PU 



1 r 



1 + (1 - 2|c 2 | 2 )(cos2^) src cos2^(^ 



(3.36) 



where (cos2^) src indicates the average of cos 26^ over the source position and energy spec- 
trum. In the usual derivations, the source term is assumed to be averaged away, yet no 
average over cos2#; is shown. However, one can get away with this in some situations that 
suppress the source term without any averaging over position or energy. 

This structure for the fully-averaged survival probability is completely general, and thus 
we interpret | C2 1 2 as the probability of hopping from one mass eigenstate to the other in the 
passage through the resonance region. Thus 



D ho P = |c 2 | = exp (-7tO) 

= exp ( -2% [ X ° dxJc 2 (x) - 2cos20„C(z) + 1 

2y/2G F EN e {ty' 2 



(3.37) 



exp —i 



,5m 
~2E 







f° dt 




h 





5m 2 



2 cos 29 v 



2V2G F EN e (t) 
5m? 



+ 1 



1/2N 



This probability characterizes the non-adiabatic nature of the evolution near the avoided 
level crossing; for purely adiabatic evolution, Ph op = 0. The limits of the integral are the 
complex turning points of Eq. (|M|), i.e., the zeros of the integrand, and are labeled such 
that Imxo > 0. The middle form is particularly convenient since then the turning points are 
located by ( = exp (±2i9 v ). This result for Phop is valid for both arbitrary mixing parameters 
and an arbitrary monotonic density profile. Since our solutions were based on the solution 
for a linear density, the form of Phop follows that for a linear density: a single exponential 
which vanishes in the adiabatic limit. 



F. Comparisons of Different Densities 



For a linear density, we must recover the linear Landau-Zener result 0,0. Equation 
( p7| ) yields 



where 



Ptp = exp (-ttQ^ 



n 



lin 



7 C 5m 2 sin 2 29 



2 AE cos20„ 



1 dN(t) 



N(t) dt 



(3.38) 



(3.39) 



as expected. 

For an exponential density, Eq. (|3.37| ) yields 



(3.40) 



where 
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(3.41) 



and 



5m 2 
~2E 



1 dN(t) 
N(t) dt 



(3.42) 



This is the leading exponential to the exact result for an exponential density [13J27]. The 
exact result is 



pexp 
hop 



exp (— 7r<5(l — cos 29 v )) — exp (— 2ti8) 
1 — exp (— 2tx5) 



(3.43) 



Equation ( |3.4U| )for the exponential density was previously obtained |E8] by connecting the 
coefficients of the coupled equations in the adiabatic basis through the complex plane |f29| . 
In Fig. (1), we compare our fully-averaged result for the survival probability in a exponential 
density (the parameters are chosen to approximate the solar density J30|) with the exact 
result. The values of the vacuum angle chosen approximate those of the two best-fit models 
for the MSW solution to the solar neutrino problem | 2"6|| . In Figs. (2) and (3), we show 
the accuracy of our approximation by comparing our source term to the exact results. The 
source term is defined as the survival probability, averaged over detector, minus the survival 
probability, averaged over both source and detector. Note that in Eq. ( |3.36| ), when Sm 2 /E is 
large, Q is large, and C2 — > 0, suppressing the source term. On the other hand, if the initial 
density is large enough, then when 5m 2 /E is small, — > tt/2, and sin 2^ — > 0, which also 
suppresses the source term. Therefore, the source term is non-zero only for intermediate 
values of 6m 2 , as illustrated in Figs. (2) and (3). 



G. Breakdown of the Mapping 

As can be seen from Fig. (1), our approximation does not hold in the extreme nonadi- 
abatic limit, where 5m 2 — ► 0. As emphasized in Ref. f5I| , the Miller-Good method only 
works well when the mapping is invertible. Given two potentials Va and Vb, the mapping is 
good only if it makes as much sense to map Va — > Vb as Vb — > Va- If this is not true, then 
the mapping is a projection, and something is lost. Invertibility may be thus be associated 
with a "sameness of topology." More precisely, when the mapping is not invertible, the 
comparison potential becomes multivalued. 

Let us consider the treatment of an exponential density. In this case, the root of the 
failure in mapping is the difference in the topology of higher-order turning points of the 
two potentials corresponding to linear and exponential densities. The turning points are 
located by ( = exp (±2i9 v ). For a linear density, there are only two turning points. For 
an exponential density, however, additional, higher-order turning points can be found by 
the transformation x — > x + 2nnx Sl where n is an integer and x s is the scale height of 
the exponential in our dimensionless units. As noted in the Appendix, we considered only 
the primary turning points, i.e., those closest to the real axis. When only the lowest-order 
turning points of the exponential density are taken into account, then the two potentials 
can be made only approximately equivalent. In principle, the way to cure this problem is to 
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use an comparison potential with the same number (infinite, if necessary) of turning points 
as the original one. In practice, this may be rather cumbersome. 

Consider how a path in the x-plane is mapped into the S'-plane. In the Appendix, we 
discuss why the locations of the primary turning points are only considered to lowest order in 
A. In particular, the turning points in the S-plane are located by S(xq) ~ So(xo) = +iy/Ti, 
S*(xo) ~ Sq(xq) = —iy/Ti. The resonance point x c is mapped to S(x c ) m 0. In the extreme 
adiabatic limit, the path from — oo to +00 along the real axis in the x-plane is mapped 
onto a path from —00 to +00 along the real axis in the S-plane. As the mixing parameters 
become more and more nonadiabatic, the path in the S'-plane makes more and more of 
an excursion into one half (upper or lower) complex plane near the resonance. At ±00, it 
returns to the real axis. In both planes, the paths run between the primary turning points. 
In the extreme nonadiabatic limit, however, the path of S(x) eventually crosses a turning 
point. There is then a topological difference between the two planes - in one case, the path 
runs between the primary turning points, and in the other, it does not. Because of how the 
turning points are anchored, this indicates that the mapping has folded the complex plane 
over, and the comparison potential is multivalued. 

The need to impose the same turning point topology between the original and comparison 
potentials restricts the applicability of Eq. ( p.37[ ) to monotonically- varying electron densities, 
i.e., those with a single MSW resonance. If there are two or more close MSW resonances, 
one cannot use a linear density to construct the comparison potential. Such situations are 
considered in Ref. p5U32f . 



IV. CONCLUDING REMARKS 

We have studied a uniform semiclassical approximation for the matter-enhanced neutrino 
oscillations for two flavors, assuming a monotonically changing but otherwise arbitrary den- 
sity profile. We obtained an analytic expression for the electron neutrino survival probability, 
unaveraged over either detector or source positions. Our result is valid for a large range of 
the mixing parameters, up to but not including the extreme nonadiabatic limit. Upon aver- 
aging over detector and source positions, we recover expressions previously obtained in the 
literature. Since our expressions are valid for arbitrary densities, they may be applied not 
only to the sun, but to all settings in which resonant neutrino conversion can occur, such as 
supernovae and the early universe. 

The method of analytic continuation utilized in Ref. [2S| for an exponential density was 



extended in Ref. |33| , where the general form of Eq. ( |3.37p for an arbitrary monotonic density 



profile was found. Results for several other analytically solvable densities are presented there. 
Our analysis not only yields an expression for the hopping probability which coincides with 
Ref. |33), but also provides the source and detector terms. 



As noted, we assumed a monotonic density profile, so this formalism is not suitable for 
studying neutrino propagation in stochastic media (e.g., with density fluctuations), as has 
recently been studied for the sun and type-II supernovae in Ref. [p2|,p4|. 
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APPENDIX: SUPERS YMMETRY-INSPIRED UNIFORM APPROXIMATION 



1. General Treatment 



Consider the Schrodinger-type equation 
d 2 ty(x) 



A 



dx 1 



A + <p 2 (x) +i\ri<p'(x) V(x)=0, (Al) 



where A is a real, positive constant, and tp is a real, monotonic function on the real axis, and 
is analytic in the complex plane. In the above and what follows, everything is dimensionless, 
and A is being used as a placeholder for h. Neither A nor (p depends on A. We will solve 
this equation in an approximation that treats A as formally small. In the physical problem 



represented by Eq. (|Al|), the variable x is real. However, for addressing the mathematical 
question of the solution of this differential equation, we consider x to be complex. We 
assume^] that A + ip 2 (x) has two zeros in the complex plane, i.e., points Xq,Xq where <p = 
±iy/A. These points are taken to be the turning points of Eq. (|A1|). We will discuss below 
why the turning points can be taken at lowest order, i.e., given as the zeros of A + ip 2 (x), 
rather than of A + (f 2 (x) + i\r]ip(x). In general, there can be more than two zeros of A + if 2 ; 
for now, we only consider the two closest to the real axis, and label them so that Imx > 0. 
The overall sign on ip is chosen to make Im<^(a;o) > 0, and r\ is taken to be ±1 as needed 
so that <p(x) has the desired sign. If presented with an equation like Eq. ( |A1| ), but with the 
opposite sign on the imaginary term, one can always conjugate it and solve as below for 
so the treatment here is general. 
We will map Eq. (|Alj) onto the comparison equation 

fi + S 2 + iXr]} U(S) = , (A2) 



A 2 d 2 U(S) 
OS 2 



where Q is a real, positive constant (and will be determined below). This equation, con- 
sidered as a function of S, also has two conjugate turning points in the complex plane. By 
"map" , we mean that we will find a change of variables S = S(x) such that the potential in 
the comparison equation is deformed into the potential in the original equation. That state- 
ment indicates how the real axis will be stretched. However, we will also have to consider 
how the complex a;-plane is mapped onto the complex S-plane. In particular, the turning 
points in the x-plane must be mapped onto the turning points in the 5-plane. The compar- 
ison equation is chosen to be one for which exact analytic solutions are known, and which 
is as similar as possible to the original equation. If we require ip(x) to be monotonic for 
real x, then imaginary term in the comparison equation may be taken as constant. Other 
than Q, this comparison equation is taken with no free parameters; such parameters can 
always be scaled away, and so are irrelevant here. That the turning-point topologies of the 
original and comparison problems be the same is critical to the method. In this case, we 



2 The case in which the zeros are on the real axis, while not relevant here, can be treated similarly 
to the rest of the Appendix; see Ref. B. 
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are mapping an as-yet unspecified real barrier onto a parabolic barrier, and the imaginary 
term onto an imaginary constant. However, we can map onto any convenient potential with 
known solutions. 

In principle, if we could find the change of variables S = S(x) exactly, then we would have 
exact solutions for ty(x) in terms of the known functions U(S(x)). In general, the solution 
for the change of variables S = S(x) would be at least as difficult as direct solution of the 
original problem. The approximation made to solve Eq. ([AI]) will be to approximate S(x) 
as a truncated power series in A. This method of uniform approximation via mapping is also 
known as the method of comparison equations (note Ref. |55|). The work here was inspired 
by the ideas of Miller and Good |IB| . For further work on the theory of their method, see 
also Ref. |36|j3ll| , and the related works in Ref. p7 |. 

In the original Miller-Good problem, the imaginary terms in the potentials above are 
not present. They treat the cases of a particle bound in a well and traveling in the presence 
of a barrier, mapping onto a parabolic well and barrier, respectively, each of which have as 
solutions the parabolic cylinder functions. In their formalism, one immediately sees that the 
WKB approximation amounts to mapping onto the free-particle potential; the mismatch in 
turning-point topologies is the origin of the failure of the WKB method (via the zeros in 
what is essentially a Jacobian, see Eq. ( |A19| ) near the turning points. With the Miller-Good 
formalism, the wave function is continuous through the turning points. 

The notation used here has some important differences from this previous work. This 
allows some difficulties and errors to be resolved. In particular, t] will be used differently 
here. We continue to take tp(x) to be real on the real axis and to be monotonic. In those 
papers, it was assumed that <p(x) monotonically increasing along the real axis would imply 
Im(p(xo) > 0. While this is suggested by the Cauchy-Riemann conditions applied to (p(x) 
on the real axis, it need not always be true. No such assumption is made here. For each 
density, one simply has to make sure that the signs of <f(x) and rj are defined so that r](p(x) 
represents the right physics and that Im<£>(xo) > 0. 

This mapping will be accomplished as 



#(x) = K(x)U(S(x)) 



(A3) 



where the form of K(x) will be chosen and S(x) will be defined by that choice. Using this 
form of ^(x) and Eq. (|A2|) , we can rewrite Eq. (|A1|) . By making the choice 



K(x) 



S'(x) 



(A4) 



and dividing through by we find 
K" 



K 



S' 



I 2 



Q + S 2 + iXr] + A + <p 2 + i\w 



So far, no approximation has been made, and the form 

1 



^(x) 



S'(x) 



-U(S(x)) 



(A5) 



(A6) 
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is a purely formal solution of Eq. (|A1|) in terms of the solutions of Eq. (|A2j). If we could find 
the change of variables S = S(x) exactly, then we would have exact solutions for in 
terms of the known functions U(S(x)). In general, the solution for the change of variables 
S = S(x) would be at least as difficult as direct solution of the original problem. To avoid 
that, we will approximate S(x). Very crudely, this procedure is a perturbation expansion 
in the shapes of the two potentials — the more they resemble each other, the more our 
approximation to the change of variables is justified, and the better our solutions will 
be. Since <p is independent of A, all of the A dependence in Eq. (|A5|) is explicit. We expand 
S(x) in powers of A: 

S(x) = S (x) + XS 1 (x) + ... . (A7) 

The power of this method is that we can obtain a good solution by keeping only the semi- 
classical terms (the lowest two orders in A). In the original Miller and Good problem [ITB 



the iXip' term was not present in the potential. Therefore, only A 2 appears in Eq. (|A5|) , and 
one can expand in A 2 instead of A, which leads to S(x) ~ Sq(x) + 0(X 2 ), making solving 
for the mapping quite simple. In our case, since A appears directly in Eq. (|A5| ), we must 
expand in A, which leads to 

S{x) ^ S {x) + XS 1 {x) , (A8) 

which makes solution of the mapping somewhat more complicated, but still much easier to 
solve than the original equation. After expansion of Eq. ( |A5| ) in A, we group by order in 
A and demand that each order vanish independently, as A is a free parameter as far as the 
mathematics are concerned. This yields the equations: 

0(A°): (A + v, 2 ) = S' 2 (n + Si) 

OiX 1 ): zV = 2(Q + Si)S' S' 1 + S' 2 (ir 1 + 2S S 1 ). K > 

While the original equation to be solved was linear, after approximation the system of 
equations to be solved is nonlinear. In particular, the equation for the 0(X 2 ) terms, which 
involves £2 is probably analytically intractable. Nevertheless, the integrations for Sq(x) 
and S\(x) can be performed, and the results are given below. In those integrations, the 
branch cut for the logarithm and square-root functions is taken along the negative real axis. 
Before solving for S (x) and Si(x), we show what will be left over. Using the relations for 
Sq(x) and Si(x), one can show 



A 



d nrtnr (x 1 



2 ^ appr 



dx 2 



A + ip 2 (x) + iXrj<p'(x) + X 2 e(x)] ^ apP r(x) = . (A10) 



In the rest of the Appendix, ty(x) always denotes the approximate wave function, and we 
drop the subscript. The degree to which the approximate solution fails to solve the exact 
differential equation is the local error, and is of the form 

e(x) =i(I) 2 "Kf)" s '° 2s! ■ 2s '° s[iir] + 2SoSi) ■ s[2(n + s2 ° ] 

+ (terms that depend on S 2 ) ■ (All) 
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The first two terms are familiar from either the WKB |38] or Miller-Good problems.^ 



The next three terms arise from the more general form of the potential considered here. 
Unfortunately, the remaining terms depend on 5*2, for which we have no analytic solution. 

The turning points of the original and comparison equations are taken to be the zeros 
of (A + <£> 2 (x)) and (Q + S 2 (x)), respectively. Since these are real for real x, the turning 
points are complex conjugates. As noted above, the turning points of original equation are 
labeled so that Imio > 0, and the sign of <p(x) is chosen so that <^(xo) = iVX. We map the 
turning points of the original equation onto the turning points of the comparison equation by 
demanding S(xo) = i\/VL. The way this correspondence is made ensures that the mapping 
does not flip the complex plane about the real axis (it is not flipped about the imaginary 
axis either, as will be shown below). These choices make it easier to avoid integration errors 
below. Note that all of the turning points are treated only at lowest order. 

The formal solution for Sq(x) can be written immediately from Eq. ( |A9|) , including the 
turning-point correspondence condition of Sq(xo) = i\/Vt: 

rx i rS {x) , 

I p (x, x ) = / dxJA + ^{x) = / dS ^Jn + SI . (A12) 

In order to evaluate Sq{x), we will first need Q, the energy of the comparison system. This 
is determined by demanding that the conjugate turning points correspond, i.e., So(xq) = 
—iyQ. When both sides of Eq. ( |A12|) are integrated between their turning points, the 
right-hand side can be done explicitly, yielding 



P dxJA + ip 2 (x) . (A13) 

Jx 



With Q determined, an implicit solution for Sq(x) can be obtained through integration of 
Eq. (|A12 ) to a general point x: 



^io) = - T + y V /n + 5 o+2 ln ^ I • (A1I) 

The fact that the solution for Sq(x) is left in this implicit form does not present any diffi- 
culties. When the asymptotic forms are used, this expression can be solved approximately. 
If the full forms of the parabolic cylinder functions are being used, then one will be taking 
a numerical approach anyway, and the solution for Sq(x) is rather easy. Using the Schwarz 
Reflection Principle and the integrals for Sq(x) and Q, one can show that Sq(x) is real for 
real x. Then I p (x, Xq) separates into real and imaginary parts as follows: 

I p (x, Xq) = Re(J p (x, Xq)) — . (A15) 



3 While these methods have the same form for the local error, the global results can be rather 
different, e.g., the transmission coefficient [16|. Note that the WKB error term is singular at the 
turning points, whereas the Miller-Good error term is bounded. 
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This will be needed to show that the exponentials in the asymptotic solution have purely 
imaginary arguments. 

To solve for Si(x), we use the first of Eq. ( |A9| ) to rewrite the second as 



17] Lf 



S[JQ + Si 



1 



We integrate from xq to x, changing variables as needed 



(iri + 2S Si) • 



(A16) 



17] 



<p(x) 



dip 



2 J v (x ) y/A + If' 1 



So(x ) Q + Q2 



(A17) 



In all three integrals, the lower limit has a positive imaginary part and the upper limit is real; 
this ensures that all of the square root functions are on the same sheet. Upon substitution 
of the limits of integration, this yields 



Si(x) 



17] 



In 



2Jn 



C2 
On 



if + VA + V 2 



+ In 



So + Jn + sz 



(A18) 



which is purely imaginary for real x, as claimed in Ref. [7|||. In those references, there were 
typographical errors (corrected here), and this was not obvious. In order to get Si right (i.e., 
pure imaginary), one has to be careful about a branch cut crossing during the integration. 
In order to facilitate that, ip is defined here with the sign noted in Eq. ( |2.5|) . 
The prefactor K(x) is given by 



K(x) 



1 



n 



Si(x) 



A + ip 2 



x 



1/4 



(A19) 



The truncation of 5" at lowest order here is consistent with our overall policy of keeping 
the two lowest orders in A, and will be justified below. This form allows us explain why the 
turning points may be taken at lowest order. Consider how the exact change of variables 
S = S(x) would work at a single point, where the local momenta in both problems could 
be taken as constant. Then the change of variables is just a scale change, and K 2 (x) is the 
ratio of the exact local momenta, the zeros of which are the exact turning points. When we 
approximate S(x) by power series in A, K 2 (x) will be the ratio of the local momenta, taken 
at the proper order, the zeros of which are the turning points, taken at the proper order. 
From the expression for K(x) above, this indicates that the turning points should be taken 
at the lowest order in A. The denominator of K(x) vanishes at the turning points. For there 
to be any hope of of ty(x) being well-behaved at the turning points, the numerator must 
vanish also. This is why we demanded the turning-point matching as above. In the WKB 
case, the numerator never vanishes, and the connection between the mismatch in turning 
point topologies and the singularity of the WKB solutions at the turning points is evident. 
When the turning points are matched properly, one can show that K(x) is well-behaved for 
all x. The method of proof is to expand Eq. ( |A12j ) near one turning point; this reveals that 
K(x) tends to a constant as the turning point is approached. The same holds for the other 
turning point as well. 
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After successively solving for Q, Sq(x), Si(x), and K(x), one can write the approximate 
wave function as 



n + s 2 (x) 

A + Lp 2 (x) 



1/4 



U(S Q {x) + XS x {x)) 



(A20) 



where U(So(x) + XSi(x)) approximately solves Eq. ( |A2|) , the comparison equation, and *ff(x) 
approximately solves Eq. ([All) , the original Schrodinger-type equation. By taking 



z(x) 



l + i 

Vx 



S(x) 



l+i 



(S (x) + XS^x)) 



(A21) 



one can show that Eq. ( |A2|) is Weber's equation ||T9|--pl|1 for the parabolic cylinder function 
D u (z), where the order v is given by 



v 



77-I -ifl/X 



The general solution may be taken to be 

*(ar) = K(x) \AD u {z{x)) + BD v {-z{x))\ . 



(A22) 



(A23) 



There are four solutions to Weber's equation, any two of which can be taken as independent; 
the two above are convenient. Given appropriate initial conditions, one can solve for ty(x) 
everywhere. The coefficients A and B will be determined below by applying the initial 
conditions to the asymptotic form of the solution. The evaluation of the parabolic cylinder 
functions is discussed in the body of the paper. 



2. Asymptotic Forms 

In many cases, one only needs the wave function as x — > ±00 (the reasons for this in 
the MSW problem are explained in the body of the paper). In this section, we develop 
the asymptotic forms of ^(x), which are easier to work with than the general form above. 
Another benefit of the asymptotic forms is that it becomes much easier to count powers of 
A, thus ensuring that we are working to a consistent order. 

Using the definitions in the previous section, one can easily show 

W^o ±0 ° (A24) 

and 

S.M-^O. (A25) 

With these crude limits, and the fact that Sq(x) is real for real x, we can determine the 
phases to the arguments of D v (z(x)) and D v (—z(x)) to be 7r/4 and — 37r/4, respectively for 
x — > +00, and vice versa for x — > —00. We use Sir/ 4 instead of 57r/4 to stay inside the 
principal branches of the square root and logarithm functions used below. In particular, one 
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must be careful when rewriting z v . The asymptotic forms of the parabolic cylinder functions 



for \z\ » \v\ are given in Ref. [20]. For Sir/A < aig(z) < 37r/4 



z 2 \ / 



D u (z) — t exp ( -- ] z v ( 1 + O ( ^ ) ) (A26) 



and for — 57r/4 < arg(z) < — 7r/4, 



2 ^ exp f ^ z-*- 1 (l + (1 ) ) . (A27) 



r(-i/) V 4 / ^ ^ : 

In the common range of validity, the difference between the two forms is negligibly small. 

These asymptotic forms make the dependence upon A in the various terms easy to see, 
and show how to keep consistent orders in A. As befits our semiclassical expansion, we only 
keep the lowest two orders in A in the exponentials of the asymptotic wave functions. Using 
the form of z given in Eq. ( |A21 ), 



exp (± Z ~Pj = exp U l - IS + 2S Sx + 0(A)] \ , (A28) 



1 + iV ^ex P (-f| + 0(A)), (A29) 



v/A 



z-^ = 1-^1 S ^ exp [+ Y ^- + 0(A) ) , (A30) 



1 + 0(1) = exp (0(A)) , (A31) 



K(x) = . = . exp (0(A)) . (A32) 
V^S) 

Now we expand the various pieces for large |x|. Using Eq. ( |A14j ), and keeping only terms 
growing or constant in |So(x)|, one can easily show 

s 2 (x) ^i P (x,x ) | (±»7T-i)n | n^/nx n ln f 2\s (x)\ \ (A33) 



2A »-±oc A 4A 4A VA/ 2A V VA 

When | So | >> ft and this expansion is valid, then \z\ >> \u\ and the D„s can be put 
in their asymptotic forms. Using Eq. (|A15 ), one can see that Sq is purely real for real x, 



as claimed. The various A terms were introduced to show where h's would appear if we 
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took these equations out of dimensionless form. The expansion for Si from Eq. ( |A18|) is 
straightforward, and yields 



SAx) 



IT] 



In 



A 



(A34) 



In this expansion, | tSo I was treated asymptotically but \tp\ was not. (At the MSW resonance 
point discussed in the body of this paper, <p(x c ) = but Sq(x c ) 7^ in general.) The 
nonadiabatic region is in general far narrower than the region over which the matter angle is 
varying appreciably. The region in which So cannot be treated asymptotically is essentially 
the nonadiabatic region, whereas (p cannot be treated asymptotically until the matter angle 
is very close to either 7r/2 or 9 V . We expand K(x) as 



K(x) 



- - 

x— >±oo \ A 



1/4 



A 



A(A + ip 2 (x)) 



1/4 



(1 ( 2\S (x) 



(A35) 



The factors of |So(ar)| remaining in these asymptotic expansions will all cancel. 
For convenience in applying the initial conditions, it is useful to write 

Ip(x, Xq) = Ip{Xi, Xq) + Ip{X, Xi) , 



(A36) 



where the first term on the right is a complex constant and the second is real and varying 
(cf. Eq. ( |A15| )). After some algebra, we find that the asymptotic expansion of the general 
solution in Eq. ( |A23j ) can be written as 



iff(x) 



x— >±oo 



A 



,1/4 



A(A + ip 2 (x)) 



±r,/2 



A 



exp ( TjI P [x,Xij 



V + ^A + v' 1 



=F*?/2 



exp ( ±-I p (x,Xi 



In the arguments of the exponentials, all terms of 0(A) or that vanish as |x| 
been dropped. The constant coefficients are given by: 



C 4 



C exp ( — — ReJ p (xj, x ) ) yA + Be 



(A37) 
oo have 

(A38) 



C- C exp [ +jReI p (xi, x )J (Ae~ iu * + B 



D + = Dexp ( +-ReI p (x i} x ) ) [Be~ 



D_ = Dexp (-jReI p (xi,Xo)^ (Ae^) 



(A39) 
(A40) 
(A41) 
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and 



= (a ) (i) bi HUa)- < A42 > 



B = -fR)U) (a) br) ««>(-«)■ < A43 > 

Note that the arguments to the x-dependent exponentials above are purely imaginary since 
I p (x,Xi) is real (that knowledge will be convenient when we take the squared modulus of 
the wave function). 
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Figure Captions 



FIG. 1. The electron neutrino survival probability vs. the mass-squared difference pa- 
rameter for two different vacuum mixing angles. The solid line is given by the method of 
this paper. The dashed line is the exact (numerical) result. The dotted line is the linear 
Landau-Zener result. In the top figure, the lines are indistinguishable. An exponential den- 
sity with parameters chosen to approximate the sun was used |30[. The region leftward of 
the lower left corner of the trough is the nonadiabatic region. 

FIG. 2. The source term (the survival probability, averaged over detector, minus the 
survival probability, averaged over both source and detector) in the electron neutrino survival 
probability vs. the mass-squared difference parameter for sin 2 29 v = 0.7. The solid line is 
given by the method of this paper. The dashed line is the exact (numerical) result. The 
density profile is as in Fig. 1. 

FIG. 3. The source term (the survival probability, averaged over detector, minus the 
survival probability, averaged over both source and detector) in the electron neutrino survival 
probability vs. the mass-squared difference parameter for sin 2 29 v = 0.01. The solid line is 
given by the method of this paper. The dashed line is the exact (numerical) result. The 
lines are indistinguishable, even when a zoom is performed in the region of rapid oscillations. 
The density profile is as in Fig. 1. 
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